clear
version 14.0
set more off
set scheme s2mono

* Specify the year of cash rent data to use in the regression
local reg_year 2013
* Specify 1 if irrigated or 0 if nonirrigated
local practice 0

capture log close
if `practice'==0 {
log using "..\logs\main results `reg_year' wheat.txt", text replace
}
if `practice'==1 {
log using "..\logs\main results irrigated `reg_year' wheat.txt", text replace
}

* run the do-files to define programs used in this code
run "..\codeSetup\soils transform - program.do"
run "main results - program.do"

local bootstrap_reps 1000

*------------------------------------------------------------------------
* Potential soils
*------------------------------------------------------------------------
local soils  slope om sar gypsum cec7 ec  ksat clay_perc silt_perc ///
			bulkDensity soc0_150 rootznemc ph_less6 ph_greater7dot5 hydgrp* text_*

local soils_nonlin  slope_* om_* sar_* gypsum_* cec7_* ec_*  ksat_* clay_perc_* silt_perc_* ///
			bulkDensity_* soc0_150_* rootznemc_* ph_less6 ph_greater7dot5 hydgrp* text_*

*------------------------------------------------------------------------
* Controls to include in the regression
*------------------------------------------------------------------------
local deficit water_deficit
local surplus water_surplus
local edd dday34C
local soils_controls soc0_150_4 bulkDensity_4 ec_2 ph_less6 ph_greater7dot5 slope_1     
*local soils_controls $SelSoils

local knots_`deficit' 2
local knots_`surplus' 2
local knots_gdd 3
local knots_`edd' 3

local main_year 0000

*------------------------------------------------------------------------
* Regression
*------------------------------------------------------------------------

use "..\dataAnalysis\rent_climate_data", clear

*******************************************************
*******************************************************
drop if winter_wheat>0.1 & winter_wheat!=.
*******************************************************
*******************************************************

if `practice'==0{
gen ln_cash_rent=ln(cash_rent_non`reg_year')
}
if `practice'==1{
gen ln_cash_rent=ln(cash_rent_irr`reg_year')
}

ricardian_program, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_s(`knots_`surplus'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			practice(`practice') reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("AVG") wild(0) 

ereturn display
matrix beta_main=e(b)

est restore main_results
	
if `reg_year'==`main_year' & `practice'==0 {
esttab main_results using "..\tables\main_results.tex", ///
			b(%9.3f) se(%9.3f) tex se label replace starlevels(* 0.10 ** 0.05) ///deficit
			r2 sfmt(%9.2f) nogaps mtitles("Coefficients") nonumbers ///
			coeflabels(	water_deficit_spline "Water Deficit" ///
						water_surplus_spline  "Water Surplus" ///
						gdd_spline1 "GDD Spline 1" gdd_spline2 "GDD Spline 2"   ///
						dday34C_spline1 "EDD Spline 1" dday34C_spline2 "EDD Spline 2" ///
						soc0_150_4 "SOC" bulkDensity_4 "Bulk Density" ///
						ec_2 "EC" ph_less6 "pH Less than 6" ph_greater7dot5 "pH Greater than 7.5" ///
						slope_1 "Log of Slope")		///
			drop(*stfips)
}


*------------------------------------------------------------------
* Graph predicted rent with average soils across water deficit
*------------------------------------------------------------------
preserve
	foreach var of varlist `surplus' gdd `edd' `soils' {
		qui summ `var', detail
		qui replace `var'=r(p50)
	}
	* Nebraska has an average cash rent very close to overall sample average
	* so assume Nebraska for the state fixed effect.
	qui replace stfips=31
	
	soils_tranform
	
	foreach var of varlist `surplus' gdd `edd' {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
			if `knots_`var''==2 {
			gen `var'_spline = `var'
			}
			if `knots_`var''==3 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
			}
			if `knots_`var''==4 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
			}
			if `knots_`var''==5 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
			}
			if `knots_`var''==6 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
			}
	}
	
	* Predict ln(rent)	
	if `practice'==0 {
	qui predict ln_rent_hat, xb 
	qui predict ln_rent_hat_se, stdp 
	}
	if `practice'==1 {
	qui predict ln_rent_hat if acresh_Irr2012!=0 & acresh_Irr2012!=., xb 
	qui predict ln_rent_hat_se if acresh_Irr2012!=0 & acresh_Irr2012!=., stdp 
	}
	* Predict rent and upper and lower 95% CI for prediction
	gen rent_hat=exp(ln_rent_hat+e(rmse)^2/2)
	gen rent_hat_low=exp(ln_rent_hat-1.96*ln_rent_hat_se+e(rmse)^2/2)
	gen rent_hat_hi=exp(ln_rent_hat+1.96*ln_rent_hat_se+e(rmse)^2/2)
	* Restrict upper bound of CI to limit the upper bound of y-axis in graph
	replace rent_hat_hi=500 if rent_hat_hi>500 & rent_hat_hi!=.
	sort `deficit'
	if `practice'==0 {
	graph twoway (rarea rent_hat_low rent_hat_hi `deficit', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_non`reg_year' `deficit', msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  `deficit', lwidth(thick) lcolor("8 81 156") lpattern(solid)), yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("Water Deficit (in)")  name(predicted_rent_deficit, replace)
	}
	if `practice'==1 {
	graph twoway (rarea rent_hat_low rent_hat_hi `deficit', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_irr`reg_year' `deficit', msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  `deficit', lwidth(thick) lcolor("8 81 156") lpattern(solid)) if ln_cash_rent!=., yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("Water Deficit (in)")  name(predicted_rent_deficit, replace)
	}
	
	if `reg_year'==`main_year' & `practice'==0 {
	graph export "..\figures\predicted_rent_deficit.pdf", replace
	
	list `deficit' rent_hat if `deficit'>14.95 & `deficit'<15.05
	list `deficit' rent_hat if `deficit'>9.95 & `deficit'<10.05
	list `deficit' rent_hat if `deficit'>4.98 & `deficit'<5.02
	
	* Knox county Missouri
	summ cash_rent_non2013 rent_hat rootznaws water_deficit ppt if stcofips==29103

	* Hancock county Illinois
	summ cash_rent_non2013 rent_hat rootznaws water_deficit ppt if stcofips==17067
	}
	
restore



*------------------------------------------------------------------
* Graph predicted rent with average soils across water surplus
*------------------------------------------------------------------
preserve
	foreach var of varlist `deficit' gdd `edd' `soils' {
		qui summ `var', detail
		qui replace `var'=r(p50)
	}
	* Nebraska has an average cash rent very close to overall sample average
	* so assume Nebraska for the state fixed effect.
	qui replace stfips=31
	
	soils_tranform
	
	foreach var of varlist `deficit' gdd `edd' {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
			if `knots_`var''==2 {
			gen `var'_spline = `var'
			}
			if `knots_`var''==3 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
			}
			if `knots_`var''==4 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
			}
			if `knots_`var''==5 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
			}
			if `knots_`var''==6 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
			}
	}
	
	* Predict ln(rent)	
	if `practice'==0 {
	qui predict ln_rent_hat, xb 
	qui predict ln_rent_hat_se, stdp 
	}
	if `practice'==1 {
	qui predict ln_rent_hat if acresh_Irr2012!=0 & acresh_Irr2012!=., xb 
	qui predict ln_rent_hat_se if acresh_Irr2012!=0 & acresh_Irr2012!=., stdp 
	}
	* Predict rent and upper and lower 95% CI for prediction
	gen rent_hat=exp(ln_rent_hat+e(rmse)^2/2)
	gen rent_hat_low=exp(ln_rent_hat-1.96*ln_rent_hat_se+e(rmse)^2/2)
	gen rent_hat_hi=exp(ln_rent_hat+1.96*ln_rent_hat_se+e(rmse)^2/2)
	* Restrict upper bound of CI to limit the upper bound of y-axis in graph
	replace rent_hat_hi=500 if rent_hat_hi>500 & rent_hat_hi!=.
	sort `surplus'
	if `practice'==0 {
	graph twoway (rarea rent_hat_low rent_hat_hi `surplus', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_non`reg_year' `surplus', msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  `surplus', lwidth(thick) lcolor("8 81 156") lpattern(solid)), yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("Water Surplus (in)")  name(predicted_rent_surplus, replace)
	}
	if `practice'==1 {
	graph twoway (rarea rent_hat_low rent_hat_hi `surplus', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_irr`reg_year' `surplus', msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  `surplus', lwidth(thick) lcolor("8 81 156") lpattern(solid)) if ln_cash_rent!=., yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("Water Surplus (in)")  name(predicted_rent_surplus, replace)
	}
	
	if `reg_year'==`main_year' & `practice'==0 {
	graph export "..\figures\predicted_rent_surplus.pdf", replace
	}
restore

*------------------------------------------------------------------
* Predicted rent across GDD
*------------------------------------------------------------------
preserve
	foreach var of varlist `deficit' `surplus' `edd' `soils' {
		qui summ `var', detail
		qui replace `var'=r(p50)
	}
	* Nebraska has an average cash rent very close to overall sample average
	* so assume Nebraska for the state fixed effect.
	qui replace stfips=31
	
	soils_tranform
	
	foreach var of varlist `deficit' `surplus' `edd' {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
			if `knots_`var''==2 {
			gen `var'_spline = `var'
			}
			if `knots_`var''==3 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
			}
			if `knots_`var''==4 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
			}
			if `knots_`var''==5 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
			}
			if `knots_`var''==6 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
			}
	}
	
	* Predict ln(rent)	
	if `practice'==0 {
	qui predict ln_rent_hat, xb 
	qui predict ln_rent_hat_se, stdp 
	}
	if `practice'==1 {
	qui predict ln_rent_hat if acresh_Irr2012!=0 & acresh_Irr2012!=., xb 
	qui predict ln_rent_hat_se if acresh_Irr2012!=0 & acresh_Irr2012!=., stdp 
	}
	* Predict rent and upper and lower 95% CI for prediction
	gen rent_hat=exp(ln_rent_hat+e(rmse)^2/2)
	gen rent_hat_low=exp(ln_rent_hat-1.96*ln_rent_hat_se+e(rmse)^2/2)
	gen rent_hat_hi=exp(ln_rent_hat+1.96*ln_rent_hat_se+e(rmse)^2/2)
	* Restrict upper bound of CI to limit the upper bound of y-axis in graph
	replace rent_hat_hi=500 if rent_hat_hi>500 & rent_hat_hi!=.
	sort gdd
	if `practice'==0 {
	graph twoway (rarea rent_hat_low rent_hat_hi gdd, fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_non`reg_year' gdd, msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  gdd, lwidth(thick) lcolor("8 81 156") lpattern(solid)), yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("GDDs")  name(predicted_rent_gdd, replace)
	}
	if `practice'==1 {
	graph twoway (rarea rent_hat_low rent_hat_hi gdd, fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_irr`reg_year' gdd, msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  gdd, lwidth(thick) lcolor("8 81 156") lpattern(solid)) if ln_cash_rent!=., yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("GDDs")  name(predicted_rent_gdd, replace)
	}
	
	if `reg_year'==`main_year' & `practice'==0 {
	graph export "..\figures\predicted_rent_gdd.pdf", replace
	}
restore


*------------------------------------------------------------------
* Predicted rent across EDD
*------------------------------------------------------------------
preserve
	foreach var of varlist `deficit' `surplus' gdd  `soils' {
		qui summ `var', detail
		qui replace `var'=r(p50)
	}
	* Nebraska has an average cash rent very close to overall sample average
	* so assume Nebraska for the state fixed effect.
	qui replace stfips=31
	
	soils_tranform
	
	foreach var of varlist `deficit' `surplus' gdd {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
			if `knots_`var''==2 {
			gen `var'_spline = `var'
			}
			if `knots_`var''==3 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
			}
			if `knots_`var''==4 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
			}
			if `knots_`var''==5 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
			}
			if `knots_`var''==6 {
			mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
			}
	}
	
	* Predict ln(rent)	
	if `practice'==0 {
	qui predict ln_rent_hat, xb 
	qui predict ln_rent_hat_se, stdp 
	}
	if `practice'==1 {
	qui predict ln_rent_hat if acresh_Irr2012!=0 & acresh_Irr2012!=., xb 
	qui predict ln_rent_hat_se if acresh_Irr2012!=0 & acresh_Irr2012!=., stdp 
	}
	* Predict rent and upper and lower 95% CI for prediction
	gen rent_hat=exp(ln_rent_hat+e(rmse)^2/2)
	gen rent_hat_low=exp(ln_rent_hat-1.96*ln_rent_hat_se+e(rmse)^2/2)
	gen rent_hat_hi=exp(ln_rent_hat+1.96*ln_rent_hat_se+e(rmse)^2/2)
	* Restrict upper bound of CI to limit the upper bound of y-axis in graph
	replace rent_hat_hi=500 if rent_hat_hi>500 & rent_hat_hi!=.
	sort `edd'
	if `practice'==0 {
	graph twoway (rarea rent_hat_low rent_hat_hi `edd', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_non`reg_year' `edd', msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  `edd', lwidth(thick) lcolor("8 81 156") lpattern(solid)), yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("EDDs")  name(predicted_rent_edd, replace)
	}
	if `practice'==1 {
	graph twoway (rarea rent_hat_low rent_hat_hi `edd', fcolor("198 219 239") fi(inten70) lcolor(white)) ///
	(scatter cash_rent_irr`reg_year' `edd', msymbol(Oh) mcolor(gs11) msize(small)) ///
	(line rent_hat  `edd', lwidth(thick) lcolor("8 81 156") lpattern(solid)) if ln_cash_rent!=., yscale(range(0 225) ) ///
	graphregion(color(white)) scale(1.5) ylabel(,nogrid angle(horizontal))  ///
	legend(off order(2 "Observed"  3 "Prediction" 1 "95% CI")) ///
	xtitle("EDDs")  name(predicted_rent_edd, replace)
	}
	
	if `reg_year'==`main_year' & `practice'==0 {
	graph export "..\figures\predicted_rent_edd.pdf", replace
	
	list `edd' rent_hat if `edd'<0.02
	list `edd' rent_hat if `edd'>2.48 & `edd'<2.52
	list `edd' rent_hat if `edd'>4.95 & `edd'<5.05
	}
	
restore

*------------------------------------
* Combine graphs
*-------------------------------------
if `reg_year'==`main_year' & `practice'==0 {
graph combine predicted_rent_deficit predicted_rent_surplus predicted_rent_gdd predicted_rent_edd, ///
		cols(2)  ycommon imargin(vsmall) graphregion(color(white)) ///
		l1title("Cash Rent ($/acre)", size(medlarge) margin(medsmall))

graph export "..\figures\predicted_rent_combined.pdf", replace	
}

if `reg_year'==`main_year' & `practice'==1 {
graph combine predicted_rent_deficit predicted_rent_surplus predicted_rent_gdd predicted_rent_edd, ///
		cols(2)  ycommon imargin(vsmall) graphregion(color(white)) ///
		l1title("Cash Rent ($/acre)", size(medlarge) margin(medsmall))

graph export "..\figures\predicted_rent_combined_irr.pdf", replace	
}

******************************************
* 20% Increase in water deficit
******************************************
preserve
	tempvar rent_hat total_rent
	estimates restore main_results
	predict ln_rent_hat, xb
	gen `rent_hat'=exp(ln_rent_hat+e(rmse)^2/2)
	drop ln_rent_hat
	gen `total_rent'=`rent_hat'*croplandNon_acres
	qui total `total_rent'
	scalar base_rent=_b[`total_rent']
	replace `deficit'=`deficit'*(1+0.2)
	damages, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_s(`knots_`surplus'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			rent_hat(`rent_hat')
	* Relative change in rent
	local perc_chg_rent=r(total_damages)/base_rent
	display `perc_chg_rent'
restore
*******************************************

******************************************
* 2 degree C increase in temp - effect on all variables
******************************************
preserve
	tempvar rent_hat total_rent
	estimates restore main_results
	predict ln_rent_hat, xb
	gen `rent_hat'=exp(ln_rent_hat+e(rmse)^2/2)
	drop ln_rent_hat
	gen `total_rent'=`rent_hat'*croplandNon_acres
	qui total `total_rent'
	scalar base_rent=_b[`total_rent']
	
	merge 1:1 stcofips using "..\temp\PRISM_county_climate83-12_GrowSeason_2C",  update replace
	keep if _merge==5
	drop _merge
	* ET and precipitation are measured in mm
	* convert ET and precipitation to inches
	local varlist `deficit' `surplus'		
	foreach var in `varlist' {
		qui replace `var'=0.0393701*`var'
	}
	replace gdd=(dday10C-dday30C)/100
	damages, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_s(`knots_`surplus'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			rent_hat(`rent_hat')
	* Relative change in rent
	local perc_chg_rent=r(total_damages)/base_rent
	display `perc_chg_rent'
restore
*******************************************

******************************************
* 2 degree C increase in temp - effect on on GDDs and EDDs
******************************************
preserve
	tempvar rent_hat total_rent
	estimates restore main_results
	predict ln_rent_hat, xb
	gen `rent_hat'=exp(ln_rent_hat+e(rmse)^2/2)
	drop ln_rent_hat
	gen `total_rent'=`rent_hat'*croplandNon_acres
	qui total `total_rent'
	scalar base_rent=_b[`total_rent']
	
	summ `edd', detail
	gen `deficit'_old=`deficit'
	gen `surplus'_old=`surplus'
	merge 1:1 stcofips using "..\temp\PRISM_county_climate83-12_GrowSeason_2C",  update replace
	keep if _merge==5
	drop _merge
	replace `deficit'=`deficit'_old
	replace `surplus'=`surplus'_old
	replace gdd=(dday10C-dday30C)/100
	summ `edd', detail
	damages, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_s(`knots_`surplus'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			rent_hat(`rent_hat')
	* Relative change in rent
	local perc_chg_rent=r(total_damages)/base_rent
	display `perc_chg_rent'
restore
*******************************************

			
simulate _b, seed(2030) reps(`bootstrap_reps') : ricardian_program, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_s(`knots_`surplus'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			practice(`practice') reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("AVG") wild(1)
if `reg_year'==`main_year' {
save "..\temp\uncertain_reg", replace 
}			

qui corr _all, cov
matrix var=r(C)

matrix colnames var = rel_total_chg rel_deficit_chg rel_surplus_chg rel_heat_chg ///
				share_deficit share_surplus share_heat 
matrix rownames var = rel_total_chg rel_deficit_chg rel_surplus_chg rel_heat_chg ///
				share_deficit share_surplus share_heat 

ereturn post beta_main var, dof(14)
ereturn display		

local total_mean_effect=_b[rel_total_chg]
local deficit_mean_effect=_b[rel_deficit_chg]
local surplus_mean_effect=_b[rel_surplus_chg]
local heat_mean_effect=_b[rel_heat_chg]




*-----------------------------------------------------------------------------------------------------
* Estimate Damages across all climate models - RCP 4.5
*-----------------------------------------------------------------------------------------------------
if `reg_year'==`main_year' & `practice'==0 {
local model_list "ACCESS bcc BNUESM CanESM2 CCSM4 CESM1BGC CNRMCM5 CSIROMk3 inmcm4 IPSLLR IPSLMR MIROC5 MIROCESM MIROCCHEM MPIESMLR MPIESMMR MRICGCM3 NorESM1M"

* dataset to store results across climate models (only climate uncertainty)
tempname memhold
tempfile results
postfile `memhold' str15 model rel_total_chg rel_deficit_chg rel_surplus_chg rel_heat_chg ///
					share_deficit share_surplus share_heat using `results', replace
* dataset to store results across boostrap reps of all climate models (regression and climate uncertainty)
tempfile bootstrap_results

foreach model of local model_list {
use "..\dataAnalysis\rent_climate_data", clear
if `practice'==0{
gen ln_cash_rent=ln(cash_rent_non`reg_year')
}
if `practice'==1{
gen ln_cash_rent=ln(cash_rent_irr`reg_year')
}
qui ricardian_program, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_s(`knots_`surplus'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			practice(`practice') reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("`model'") wild(0)

ereturn display
matrix beta=e(b)
post `memhold' ("`model'") (beta[1,1]) (beta[1,2])	(beta[1,3])	(beta[1,4])	(beta[1,5]) (beta[1,6])	(beta[1,7])

simulate _b, seed(2030) reps(`bootstrap_reps') : ricardian_program, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_`deficit'') knots_s(`knots_`surplus'') knots_g(`knots_gdd') knots_e(`knots_`edd'') ///
			practice(`practice') reg_year(`reg_year') soils_controls(`soils_controls') scenario("rcp45") model("`model'") wild(1)

capture append using `bootstrap_results'
save `bootstrap_results', replace
			
qui corr _all, cov
matrix var=r(C)

matrix colnames var = rel_total_chg rel_deficit_chg rel_surplus_chg rel_heat_chg ///
				share_deficit share_surplus share_heat
matrix rownames var = rel_total_chg rel_deficit_chg rel_surplus_chg rel_heat_chg ///
				share_deficit share_surplus share_heat 

ereturn post beta var, dof(14)
ereturn display		

}

postclose `memhold'
use `results', clear
summ, detail
save "..\temp\uncertain_climate", replace

use `bootstrap_results', clear
save "..\temp\uncertain_climate_reg", replace

****************************
* Statistics in the paper
****************************
use "..\temp\uncertain_climate", clear
sort rel_total_chg
list model rel_total_chg

sort rel_deficit_chg
list model rel_deficit_chg

sort rel_heat_chg
list model rel_heat_chg
} // end of if `reg_year'==`main_year'

log close



if `reg_year'==`main_year' & `practice'==0 {
*-------------------------------------------------------------------
* Create graphs of uncertainty in climate change damages
*-------------------------------------------------------------------
set more off
use "..\temp\uncertain_reg", clear
foreach j in total deficit surplus heat {
	summ _b_rel_`j'_chg, detail
	foreach q in 50 75 25 95 5 {
		scalar `j'_reg_p`q'=r(p`q')
	}
}

use "..\temp\uncertain_climate", clear
foreach j in total deficit surplus heat {
	summ rel_`j'_chg, detail
	foreach q in 50 75 25 95 5 {
		scalar `j'_climate_p`q'=r(p`q')
	}
}


use "..\temp\uncertain_climate_reg", clear
foreach j in total deficit surplus heat {
	summ _b_rel_`j'_chg, detail
	foreach q in 50 75 25 95 5 {
		scalar `j'_both_p`q'=r(p`q')
	}
}


clear
set obs 3
gen uncertainty=_n
foreach q in 50 75 25 95 5  {
	gen total_p`q'=. 
	gen deficit_p`q'=. 
	gen surplus_p`q'=. 
	gen heat_p`q'=. 
}


foreach j in total deficit surplus heat {

* Regression uncertainty
foreach q in 50 75 25 95 5  {
replace `j'_p`q'=`j'_reg_p`q' if uncertainty==1
}


* Climate uncertainty
foreach q in 50 75 25 95 5  {
replace `j'_p`q'=`j'_climate_p`q' if uncertainty==2
}

* Climate & Regression uncertainty
foreach q in 50 75 25 95 5  {
replace `j'_p`q'=`j'_both_p`q' if uncertainty==3
}

* Format the mean effects to display on the graphs
local `j'_mean_effect : di %9.2f ``j'_mean_effect'
}

/*
* Relative change in total rent
twoway rspike total_p75 total_p95 uncertainty, pstyle(p1) || ///
rspike total_p25 total_p5 uncertainty, pstyle(p1) || ///
rbar total_p50 total_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar total_p50 total_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
xtitle("Source of Uncertainty", size(medlarge) margin(medsmall)) ///
graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
text(-0.25 0.95 "`total_mean_effect'", place(center)) name(rel_chg_total)

graph export "..\figures\rel_chg_total.pdf", replace

* Relative change due to deficit
twoway rspike deficit_p75 deficit_p95 uncertainty, pstyle(p1) || ///
rspike deficit_p25 deficit_p5 uncertainty, pstyle(p1) || ///
rbar deficit_p50 deficit_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar deficit_p50 deficit_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
xtitle("Source of Uncertainty", size(medlarge) margin(medsmall)) ///
graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
text(-0.05 0.95 "`deficit_mean_effect'", place(center)) name(rel_chg_deficit)

graph export "..\figures\rel_chg_deficit.pdf", replace

* Relative change due to surplus
twoway rspike surplus_p75 surplus_p95 uncertainty, pstyle(p1) || ///
rspike surplus_p25 surplus_p5 uncertainty, pstyle(p1) || ///
rbar surplus_p50 surplus_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar surplus_p50 surplus_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
xtitle("Source of Uncertainty", size(medlarge) margin(medsmall)) ///
graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
text(0.02 0.95 "`surplus_mean_effect'", place(center)) name(rel_chg_surplus)

graph export "..\figures\rel_chg_surplus.pdf", replace

* Relative change due to heat
twoway rspike heat_p75 heat_p95 uncertainty, pstyle(p1) || ///
rspike heat_p25 heat_p5 uncertainty, pstyle(p1) || ///
rbar heat_p50 heat_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar heat_p50 heat_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
xtitle("Source of Uncertainty", size(medlarge) margin(medsmall)) ///
graphregion(color(white)) scale(1.5) ylabel(,nogrid) ///
text(-0.15 0.95 "`heat_mean_effect'", place(center)) name(rel_chg_heat)

graph export "..\figures\rel_chg_heat.pdf", replace
*/


* Relative change in total rent
twoway rspike total_p75 total_p95 uncertainty, pstyle(p1) || ///
rspike total_p25 total_p5 uncertainty, pstyle(p1) || ///
rbar total_p50 total_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar total_p50 total_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
title("Total") ///
graphregion(color(white)) scale(1.25) ylabel(,nogrid) ///
text(-0.2 0.95 "`total_mean_effect'", place(center)) name(rel_chg_total, replace)

graph export "..\figures\rel_chg_total.pdf", replace

* Relative change due to deficit
twoway rspike deficit_p75 deficit_p95 uncertainty, pstyle(p1) || ///
rspike deficit_p25 deficit_p5 uncertainty, pstyle(p1) || ///
rbar deficit_p50 deficit_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar deficit_p50 deficit_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
title("Water Deficit") ///
graphregion(color(white)) scale(1.25) ylabel(,nogrid) ///
text(0.0 0.95 "`deficit_mean_effect'", place(center)) name(rel_chg_deficit, replace)

graph export "..\figures\rel_chg_deficit.pdf", replace

* Relative change due to surplus
twoway rspike surplus_p75 surplus_p95 uncertainty, pstyle(p1) || ///
rspike surplus_p25 surplus_p5 uncertainty, pstyle(p1) || ///
rbar surplus_p50 surplus_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar surplus_p50 surplus_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
title("Water Surplus") ///
graphregion(color(white)) scale(1.25) ylabel(,nogrid) ///
text(0.12 0.95 "`surplus_mean_effect'", place(center)) name(rel_chg_surplus, replace)

graph export "..\figures\rel_chg_surplus.pdf", replace

* Relative change due to heat
twoway rspike heat_p75 heat_p95 uncertainty, pstyle(p1) || ///
rspike heat_p25 heat_p5 uncertainty, pstyle(p1) || ///
rbar heat_p50 heat_p75 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35) || ///
rbar heat_p50 heat_p25 uncertainty, pstyle(p1) bfc(gs15) blc(gs8) barw(0.35)  ///
legend(off) xla(1 "Regression" 2 "Climate" 3 "Both", noticks) ///
title("Heat Stress") ///
graphregion(color(white)) scale(1.25) ylabel(,nogrid) ///
text(-0.05 0.95 "`heat_mean_effect'", place(center)) name(rel_chg_heat, replace)

graph export "..\figures\rel_chg_heat.pdf", replace

graph combine rel_chg_total rel_chg_deficit rel_chg_surplus rel_chg_heat, ///
		cols(2)  xcommon ycommon imargin(vsmall) graphregion(color(white)) ///
		b1title("Source of Uncertainty", size(medlarge) margin(medsmall)) ///
		l1title("Relative Change in Rent", size(medlarge) margin(medsmall))

graph export "..\figures\rel_chg_combined.pdf", replace	
		


} // end of if `reg_year'==`main_year'

